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^ : ABSTRACT 

O ! We consider one-loop corrections (non-linear corrections beyond leading or- 

Q ! der) to the bispectrum and skewness of cosmological density fluctuations induced 

^ ! by gravitational evolution, focusing on the case of Gaussian initial conditions and 

CN ! scale-free initial power spectra, P{k) oc k"'. As has been established by compar- 

^ ! ison with numerical simulations, tree-level (leading order) perturbation theory 

^ ! describes these quantities at the largest scales. One-loop perturbation theory 

O i provides a tool to probe the transition to the non-linear regime on smaller scales. 

^ I In this work, we find that, as a function of spectral index n, the one-loop bispec- 

T— I I trum follows a pattern analogous to that of the one-loop power spectrum, which 

^ I shows a change in behavior at a "critical index" Uc ~ —1-4, where non-linear 

^ I corrections vanish. For the bispectrum, for n < one-loop corrections increase 

the configuration dependence of the leading order contribution; for n > Uc, one- 
O I loop corrections tend to cancel the configuration dependence of the tree-level 

I bispectrum, in agreement with known results from n = —1 numerical simula- 

I tions. A similar situation is shown to hold for the Zel'dovich approximation, 

I where ric ~ —1.75. Using dimensional regularization, we obtain explicit analytic 

^ I expressions for the one-loop bispectrum for n = —2 initial power spectra, for 

I both the exact dynamics of gravitational instability and the Zel'dovich approx- 

imation. We also compute the skewness factor, including local averaging of the 
density field, for n = —2: SslR) = 4.02 + 3.83 (Tq{R) for gaussian smoothing and 
Ss{R) = 3.86 + 3.18 a'^^{R) for top-hat smoothing, where cr^{R) is the variance 
of the density field fluctuations smoothed over a window of radius R. Comparison 
with fully non-linear numerical simulations implies that, for n < —1, one- loop 
perturbation theory can extend our understanding of nonlinear clustering down 
to scales where the transition to the stable clustering regime begins. 
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1. Introduction 

There is growing evidence that the large-scale structure of the Universe grew via grav- 
itational instability from small primordial fluctuations in the matter density. For realistic 
models of structure formation, the initial spectrum of perturbations is such that at large 
scales, fluctuations are small and reflect the primordial spectrum. The variance of density 
fluctuations, (7^(i?), is a decreasing function of scale R. At small scales, (7^(i?) is large 
enough that non-linear effects become important. There are therefore two limiting regimes 
characterized by the value of cr^(i?): the linear regim,e at large scales, where cr^{R) <S 1, 
and the non-linear regim,e at small scales, where cr^(i?) ^ 1. The boundary between these 
two regimes defines a length scale, the correlation length i?o, where o"^(i?o) = 1- Because of 
gravitational instability, Rq grows with time and therefore a given scale eventually becomes 
non-linear under time evolution. 

At early epochs, the growth of density perturbations can be described by linear per- 
turbation theory, provided that the linear power spectrum P{k) falls off less steeply than 
/c^ for small k (Zel'dovich 1965, Peebles 1974, Peebles & Groth 1976). In the linear regime, 
perturbation Fourier modes evolve independently of one another, conserving the statistical 
properties of the primordial fluctuations. In particular, if the primordial fluctuations are 
Gaussian random fields, they remain Gaussian in linear theory. In this case, the statisti- 
cal properties of the density and velocity fields are completely determined by the two-point 
correlation function or the power spectrum. 

When the fiuctuations become non-linear, coupling between different Fourier modes 
becomes important, inducing non-trivial correlations that modify the statistical properties of 
the cosmological fields. For Gaussian initial conditions, this causes the appearance of higher- 
order reduced correlations, which constitute independent statistics that can be measured in 
observational data and numerical simulations, even when the departure from the linear 
regime is small. 

Non-linear cosmological perturbation theory provides a theoretical framework for the 
calculation of the induced higher-order correlation functions in the weakly non-linear regime, 
defined by scales R such that o"(i?) < 1. At large scales, leading order (tree- level) pertur- 
bation theory gives the first non-vanishing contribution, and has been used to understand 
the generation of higher order correlations in gravitational instability. Comparison with fully 
non-linear numerical simulations has shown this approach to be very successful (Juszkiewicz, 
Bouchet & Colombi 1993, Bernardeau 1994b, Lokas et al. 1995, Gaztanaga & Baugh 1995, 
Baugh, Gaztanaga & Efstathiou 1995). 

As one approaches smaller scales, however, next to leading order (loop) corrections 
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to the tree- level results are expected to become important. The question then arises of 
whether our understanding of non-linear clustering can be extended from the largest scales 
into the transition region to the non-linear regime. This motivates us to consider one-loop 
cosmological perturbation theory. In previous work (Scoccimarro & Frieman 1996b, hereafter 
SF2), we showed that for scale-free initial conditions, P{k) oc /c", without too much small- 
scale power (spectral index n < —1), one can understand the evolution of the power spectrum 
down to scales where it begins to go over to the strongly non-linear stable clustering regime. 
Therefore, it is interesting to consider one-loop corrections to the higher order correlation 
functions as well, to see if one can gain similar understanding of non-linear clustering on 
intermediate scales. 

In this work we concentrate on one- loop corrections to the three-point function of density 
perturbations in Fourier space, also known as the bispectrum, and its one-point counterpart, 
the skewness. These are interesting quantities for several reasons. On the theoretical side, 
the bispectrum is the lowest order correlation function which, for Gaussian initial conditions, 
vanishes in the linear regime; its structure therefore reflects truly non-linear properties of the 
matter distribution. Furthermore, as the lowest order correlation function which depends on 
the vector character of its arguments, it gives direct physical information on the anisotropic 
structures and flows generated by gravitational instability. Observationally, the configuration 
dependence of the tree-level bispectrum has been put forward as a promising statistic to 
study the important but poorly understood issue of bias (Fry 1994), i.e., the degree to which 
luminous objects in the universe such as galaxies are fair tracers of the underlying density 
field. It is therefore important to see how further non-linear effects (which are inevitably 
present in observational studies) alter this configuration dependence, to check whether one 
can still disentangle nonlinear evolution from bias. 

We focus on Gaussian initial conditions and scale- free initial power spectra, P{k) ex k'^. 
In addition to mathematical simplicity, primordial Gaussian fluctuations have a broad phys- 
ical motivation and are predicted by the simplest inflationary models. Although the linear 
power spectrum for the Universe is not scale-free (on both observational and theoretical 
grounds), scale- free spectra are very useful approximations over limited ranges of wavenum- 
bcr k. They also have the advantage of yielding analytic closed form results and giving rise 
to self-similar evolution of the statistical properties of cosmological fields for spatially flat 
universes (Davis & Peebles 1977, Peebles 1980). In particular, self-similarity is a powerful 
aid towards a physical understanding of non-linear clustering and in many realistic models 
of structure formation we expect approximate self-similar evolution over a restricted range 
of length and time-scales (Efstathiou et al. 1988). 

While the agreement between tree-level perturbation theory and numerical simulations 
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in the weakly non-linear regime is well established, it was not until recent years that N-body 
simulations have been able to reliably follow the transition of higher order statistics into 
the non-linear regime. In this regard, for scale-free initial power spectra, the bispectrum in 
numerical simulations has been shown by Fry, Melott & Shandarin (1993,1995) to depart 
from the tree-level perturbative results at scales comparable to the correlation length, as 
expected if next to leading order corrections are present. Similarly, for the skewness of 
the density field, and for higher-order cumulants as well, deviations from the leading order 
calculations have been reported in the literature (Bouchet & Hernquist 1992, Lucchin et 
al. 1994, Juszkiewicz et al. 1995, Hivon et al. 1995, Colombi, Bouchet & Hernquist 1996). 
We therefore consider it appropriate to extend the leading order calculations to one- loop, 
in order to understand better the limitations of the tree-level results and see the extent to 
which one can improve the agreement of perturbation theory with fully non-linear numerical 
simulations. 

This paper is organized as follows. In Section || we discuss the cosmological fluid equa- 
tions of motion and their solution within the framework of perturbation theory. For re- 
cent reviews of perturbation theory see Bernardeau (1996), Bouchet (1996), Juszkiewicz & 
Bouchet (1996); approximation methods in gravitational clustering are reviewed by Sahni & 
Coles (1996). Section |^ reviews the diagrammatic approach to perturbation theory and the 
self- similarity properties of statistical quantities derived from it. The main results of this 
work are presented in Section ^, where we consider for the first time one-loop corrections 
to the bispectrum and skewness including smoothing effects. We compare the latter with 
results from numerical simulations; a similar comparison for the bispectrum will be presented 
elsewhere (Scoccimarro et al 1996). Section |^ contains our conclusions. Auxiliary material 
is consider in the Appendices. 



The equations of motion relevant to gravitational instability describe conservation of 
mass and momentum and the Poisson equation for a self-gravitating perfect fluid with zero 
pressure in a homogeneous and isotropic universe (Peebles 1980): 



2. Dynamics and Perturbation Theory 



2.1. Equations of Motion 



d6{^, t) 
dr 



+ V-{[l + 5(x,r)]v(x,r)} = 0, 



(1) 



dv{x, t) 



+ n{T) v(x, r) + [v(x, r) ■ V]v(x, r) 



V$(x,r), 



(2) 



dr 
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V2$(x,r) = ^^]7^2(^)5(x,r) (3) 

Here, x denotes comoving spatial coordinates, r = J dt/ a is the conformal time, a(r) is 
the cosmic scale factor, the density contrast 5(x, r) = p{-x,T)/p — 1, with p(r) the mean 
density of matter, v = d'^/dr represents the velocity field fluctuations about the Hubble 
flow, Ti = dlna/dr = Ha is the conformal expansion rate, $ is the gravitational potential 
due to the density fluctuations, and the density parameter Q = p/ pc = SuGpa"^ /3H'^. We 
take the velocity field to be irrotational, so it can be completely described by its divergence 
^ = V ■ V. We will refer to Eqs. (|I|)-(0) as the "exact dynamics" (ED), to make a dis- 
tinction with the modified dynamics introduced by non-linear approximations such as the 
Zel'dovich approximation and the Local Lagrangian approximation to be discussed later (see 
Appendices ^ and 0). Equations (|l])-(^ hold in an arbitrary homogeneous and isotropic 
background Universe which evolves according to the Friedmann equations; henceforth, for 
simplicity we assume an Einstein-de Sitter background, f2 = 1, with vanishing cosmological 
constant, for which a cx and 3f27i^/2 = 

Taking the divergence of Equation (|^) and Fourier transforming the resulting equations 
of motion we get: 



dSjk, r) 

dr 



e{^, r) = -J d'h J rf'A;25D(k - ki - k2)a(k, ki)^(ki, r)5(k2, r), (4) 



- J d'h J d'hSnik - ki - k2)/3(k, ki, k,)9{k^, T)9{k,, r), (5) 
{6d denotes the three-dimensional Dirac delta distribution), where the functions 



a(k, ki) = — ^, /3(k, ki, k2) = (6) 



encode the non-linearity of the evolution (mode coupling) and come from the non-linear 
terms in the continuity equation (|I|) and the Euler equation (^ respectively. 



2.2. Perturbation Theory Solutions 



We focus on a statistical description of cosmological perturbations: we are interested in 
correlation functions of the fields 5(k, r) and 6'(k, r) (i.e., the ensemble average of products 
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of these fields). Ensemble averaging effectively introduces a new parameter into the problem, 
the variance of the density fluctuations cr^ =< 5^ >, which controls the transition from the 
linear (cr^ ^ 1) to the non-linear regime (cr^ ^ 1). We consider perturbations about 
the linear solution, effectively treating the variance of the linear fluctuations as a small 
parameter. In this case, Eqs. (H)-® can be formally solved via a perturbative expansion. 



5(k, t) = Y: a^{r)Sn{k), e{k, t) = H{t) ^ a"(r)^„(k), (7) 

n=l 71=1 

where only the fastest growing mode at each order is taken into account. At small a, 
the series are dominated by their first terms, and since ^i(k) = — 5i(k) from the continuity 
equation, (5i(k) completely characterizes the hnear fluctuations. The equations of motion (|[)- 
(|) determine (5n(k) and 6'„(k) in terms of the linear fluctuations, 

5„(k) = j £qi ... j (fqJoik - qi - ... - q„)F^"nqi, • • • , qn)5i(qi) . . . 5i(q„), (8) 
e„(k) = - j Sq, ... 1 SqjD{k - qi - ... - q„)GW(qi, . . . , q,)5i(qi) . . . 5i(q„), (9) 

where F^^ and G^*^ are symmetric homogeneous functions with degree zero of the wave 
vectors {qi, . . . , q„}. They are constructed from the fundamental mode coupling functions 
a(k, ki) and /3(k, ki, k2) according to the recursion relations {n > 2, see Goroff et al. (1986) 
or Jain & Bertschinger (1994) for a derivation): 



i:^ / \ \^ ^m(qi' ■ ■ ■ ; q^ 

F„(qi,...,q„) = 2^ 

m 



\ (2n + 3)(n- 1) 
+ 2/5(k,ki,k2)G 

n—m\S\.m+l^ ■ ■ ■ ■> ^.n) 



(2n + l)a(k, ki)Fn-m{qm+l, • • • , qn) 



(10) 



/ \ G^7Ti(qi) • • • 5 qi 

G„(qi,...,q„) = 2^ 

m 



[ (2n + 3)(n- 1) 
+ 2n/5(k,ki,k2)G' 

n—m\^.m+l: • • • ; q7ij 



3a (k, ki)F„,_^(q^+i, . . . , q„) 



(where ki = qi + . . . + q^, k2 = qm+i + . . . + q„, k = ki + k2, and ^1 = 6*1 = 1) and the 
symmetrization procedure: 



qi,..,qr, 



1 



q7r(l), • • • , q7r(7i)J 



(12) 
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Gl^\qi,...,qn) = -rI^G'n(q^(i),...,q^(„)), (13) 
n. ^ 

where the sum is taken over all the permutations tt of the set {1, . . . , n}. 

3. Statistics and Diagrammatics 

3.1. Diagrammatic Expansion of Statistical Quantities 

The starting point for a statistical description of fluctuations in cosmology is the "Fair 
Sample Hypothesis" (Peebles 1980, Bertschinger 1992). This asserts that fluctuations can 
be described by statistically homogeneous and isotropic random fields (so that our Universe 
is a random realization from a statistical ensemble) and that within the accessible part of 
the Universe there are many independent samples that can be considered to approximate a 
statistical ensemble, so that spatial averages are equivalent to ensemble averages ("ergodic- 
ity"). In this work we focus on the non-linear evolution of the three-point cumulant of the 
density field, the bispectrum i?(ki,k2,r), and its 1-point counterpart, the skewness factor 
S3{R,t). These are defined respectively by: 

(5(ki, r)5(k2, r)5(k3, r))^ = 5^(ki + ka + k,) B{k,, k,, r), (14) 

and 

S3{R,r) = ,,l . I S(ki,k2,T) W{k^R)W{k2R)W{\k^ + k2\R) d^kA, (15) 

<T [K, T) J 

where the angle brackets denote ensemble averaging , the subscript "c" stands for the con- 
nected contribution (see below), and (7^(i?, r) is the variance of the density field fluctuations: 

a\R, t) = J P{k, t) W\kR) d^k = (6\R, r)). (16) 
Here the power spectrum P{k, r) is defined by 



and therefore 



{S{k, T)~6{k', r))^ = 6D{k + k')P{k, r). 



(17) 



8 



S3{R,t) 



(18) 



Here W{kR) is the Fourier transform of the window function, which we take to be either a 
top-hat (TH) or a Gaussian (G), 



It is convenient to define the hierarchical amphtude Q as follows (Fry & Seldner 1982, Fry 
1984): 



which has the desirable property that it is scale and time independent to lowest order (tree- 
level) in non-linear perturbation theory. In a pure hierarchical model, Q would be a fixed 
constant, independent of configuration and of the power spectrum P{k,T) as well. 

We are interested in calculating the non-linear evolution of these statistical quantities 
from Gaussian initial conditions in the weakly non-linear regime, cr{R) ^1. A systematic 
framework for calculating correlations of cosmological fields in perturbation theory has been 
formulated using diagrammatic techniques (GorofT ct al. 1986, Wise 1988, Scoccimarro & 
Frieman 1996 (SFl) , SF2). In this approach, contributions to p-point cumulants of the 
density field come from connected diagrams with p external (solid) lines and r — p — l,p, . . . 
internal (dashed) lines. The perturbation expansion leads to a collection of diagrams at 
each order, the leading order being tree-diagrams, the next to leading order 1-loop diagrams 
and so on. In each diagram, external lines represent the spectral components of the fields 
we arc interested in (e.g., 5(k, r)). Each internal line is labeled by a wave- vector that is 
integrated over, and represents a linear power spectrum Pii(g, r). Vertices of order n (i.e., 
where n internal lines join) represent an n*^ order perturbative solution 6n, and momentum 
conservation is imposed at each vertex. 

We can write the loop expansion for the power spectrum up to one-loop corrections as 




(19) 



Wg{u) = exp(-'wV2). 



(20) 



Q = 




(21) 



(22) 
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where the superscript (n) denotes an n-loop contribution, the tree-level (0-loop) contribution 
is just the hnear spectrum, 



(23) 



with a^(r)(5i(k)5i(k'))c = S^ik + k')Pii(A;,r), and the 1-loop contribution consists of two 
terms, 



with (see Fig. |l]): 



p(i)(A;,r) = P22(fc,r) + Pi3(fc,r) 



(24) 



P22{k,T) = 2J[Fi'\k-q,q)]^Pu{\k-q\,T)P^,{q,r)d\ 
Pi3{k,T) ^ gJ Fi'\k,q,-q)P^,{k,r)P^^iq,T)d\ 



(25) 
(26) 



Here Pij denotes the amplitude given by the above rules for a connected diagram representing 
the contribution from {6i6j)c to the power spectrum. We have assumed Gaussian initial 
conditions, for which Pij vanishes if i + j is odd. 

For the smoothed variance we write 



(t2(P) =a|(P)(l + s« a2(P) 



(27) 



+ 



k — q 



' ' k 



11 



{P. 



22 



(A3) 



Fig. 1. — Diagrams for the power spectrum up to one-loop. See Eqs. ( pSD and (^) for 
diagram amplitudes. 



where (jj{R) denotes the variance in hnear theory (given by ([T6|) with P = Pn); the dimen- 
sionless 1-loop amphtude is 

s^'\R) = ^ / P^'\k,T) W\kR) d'k 

To characterize the degree of non-hnear evolution when including one-loop corrections 
to the power spectrum and bispectrum, it is convenient to define a physical scale from the 
linear power spectrum, the correlation length Rq, as the scale where the smoothed linear 
variance is unity, 

aj{Ro) = J d'k Pn{k,T) W\kRo) = 1. (29) 
The loop expansion for the bispectrum reads: 

5(ki, k2, r) = (ki, k2, r) + (ki, k2, r) + . . . , (30) 

where the tree-level part is given by a single diagram in second order perturbation theory 
(see Fig. ||) plus its permutations over external momenta (recall that ki + k2 -|- ka = 0): 

S(°)(ki,k2,r) = 2Pn{kur)Pn{k2,T)Fi'\kuk2) + 2Pn{k2,r)Pn{ks,r) 

xFi^)(k2, ks) + 2Pn(fc3, r)Pu{ki, T)F^'\k^, k,). (31) 

The one-loop contribution consists of four distinct diagrams involving up to fourth-order 
solutions: 

5(i)(ki, k2, r) = 5222(ki, k2, r) + 53^2i(ki, k2, r) + B^,{ki, k^, r) + 54ii(ki, k^, r), (32) 
where: 

8 J rfVii(g,T)Fi^Vq,q + ki)Pii(|q + ki|,r)Fi^Vq-ki,q-k2) 

xPn(|q-k2|,r)Fi^)(k2-q,q), (33) 
QPn{k3,T) J (-q,q-k2,-k3)Pii(|q-k2|,r) 



(28) 



222 



-°321 — 
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xF2*-'(q, k2 — q) + permutations, (34) 

B^, ^ 6Pn(^2,r)Pn(fc3,r)Fi^nk2,k3)|rfVii(g,r)Fi^nk3,q,-q) 

+permutations, (35) 

B,u = 12Pii(A;2, T)Pu{k3, r) J d\Pn{q, r)Fi^)(q, -q, -ka, -kg) 

+permutations. (36) 

For the hierarchical amphtude Q (see Eq. (^)), the loop expansion yields: 

^ ^ i?W(ki,k2,r)+pW(ki,k2,r) + ... 

S(0)(ki,k2,r) + S«(ki,k2,r) + ...' ^-"^ 

where: 



E(°)(ki, k2, r) = Pn(A;i, r)Pn(A;2, r) + Pn(A;2, r)Pn(A;3, r) + Pn(fc3, r)Pn(A;i, r), (38) 

and: 

S(^)(ki, k2, r) = P(°H^1> r)P^^\k2, r) + permutations. (39) 
For large scales, it is possible to expand Q = Q^^^ + Q^^^ + . . ., which gives: 

n(o) = ^^°Hki,k2,r) 

-S(0)(ki,k2,r)' ^^^^ 

^(1) ^ P(^(ki,k2,r)-gWsW(ki,k2,r) ^ ~ _ S(^)(ki, k2, r) 
- S(o)(ki,k2,r) S(0)(k„k2,r)- 

Note that Q*-^-* depends on the normalization of the linear power spectrum, and its amplitude 
increases with time evolution. On the other hand, from Equations (^) and ( ^O]) it follows 
that Q^^^ is independent of time and normalization (Fry 1984). Furthermore, for scale-free 
initial conditions, Pu{k) oc fc", Q^^^ is also independent of overall scale. For the particular 
case of equilateral configurations {ki = k2 = k^ and ki ■ kj = —0.5 for all pairs), Q^^^ is 
independent of spectral index as well, Q^^q = 4/7. In general, for scale- free initial power 
spectra, Q^^^ depends on configuration shape through, e.g., the ratio ki/k2 and the angle 
9 defined by ki ■ k2 = cos 6*. Note that we also defined Q^^^ in Eq. (^T]) which denotes the 



(^211) 

Fig. 2. — Tree-level diagram for the bispectrum. This diagram plus its 2 permutations over 
external momenta generates the tree-level bispectrum. The corresponding amplitudes are 
given by Eq. 



q + k{ 



q - k2i + 

k3 



+ 



(-B222) 



(-S321) 



(S411) 



Fig. 3. — One-loop diagrams for the bispectrum. The corresponding amplitudes are given 
in Eqs. (|33|) through (|36D . 
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Fig. 4. — The tree-level hierarchical amplitude Q^^^ for triangle configurations given by 
= 2 as a function of the angle 9 (fci ■ k2 = cos 9). The different curves correspond to 
spectral indices n = —2, —1.5, —1, —0.5, (from top to bottom). See also Fry (1994) 

one-loop correction to the bispectrum normalized by the tree-level quantity S^^-'; this will be 
useful in order to assess the behavior of the one-loop bispectrum with spectral index. 

Figure ^ shows Q^^^ for the triangle configuration given by ki/k2 = 2 as a function of 9 
for different spectral indices. The configuration dependence of Q^^^ is remarkably insensitive 
to other cosmological parameters, such as the density parameter Q and the cosmological 
constant (Fry 1994) (see also Hivon et al. (1995)). In fact, since bias between the galaxies 
and the underlying density field is known to change this configuration dependence (Fry 
& Gaztanaga 1993), measurements of the hierarchical amplitude Q in galaxy surveys could 
provide a measure of bias which is insensitive to other poorly known cosmological parameters 
(Fry 1994), unlike the usual determination from peculiar velocities which has a degeneracy 
with the density parameter Q. 

The configuration dependence of comes from the second order perturbation theory 
kernel ^2'^^ (see Eqs. (|40|) and (|3l|) ) and can be understood in physical terms as follows. 
From the recursion relations given in Eq. (0), we can write: 
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Fi'\k,M) = ^ f«(k,ki) + «(k,k2)l + ^ /5(k,ki,k2), (42) 



2 
7 



where k = ki + k2 (see Eq. (|^) for definitions of the mode- coupling functions a and /5). The 
terms in square brackets contribute a constant term, independent of configuration, coming 
from the 6 x S term in the equations of motion, plus terms which depend on configuration 
and describe gradients of the density field in the direction of the flow (i.e., the term v ■ VS 
in the continuity equation). Similarly, the last term in Eq. ( |i2| ) contributes configuration 
dependent terms which come from gradients of the velocity divergence in the direction of the 
flow (due to the term (v- V)v in Euler's equation). Therefore, the configuration dependence 
of the bispectrum reflects the anisotropy of structures and flows generated by gravitational 
instability. The enhancement of correlations for coUinear wavevectors {6 = 0, vr) in Figure ^ 
reflects the fact that gravitational instability generates density and velocity divergence gra- 
dients which are mostly parallel to the flow. Upon ensemble averaging, which by ergodicity 
corresponds to weighting configurations by their number frequency, this leads to a predomi- 
nance of correlations in nearly collinear configurations. The dependence on the spectrum is 
also easy to understand: models with more large-scale power (smaller spectral indices n) give 
rise to anisotropic structures and flows with larger coherence length, which upon ensemble 
averaging leads to a more anisotropic bispectrum. We will see in the next Section that this 
physical picture provides some insight into the behavior of one-loop corrections. 

The loop expansion for the skewness factor gives (SFl): 

_ Sj'^ + si%\R) + . . . 
^'^^^ = l + 2.«a2(i?) + ... ' 

where: 

5f (^) = J d'h£hB^''\k,,k2) WihR)W{k2R)W{\k, + k2\R), (44) 

Si'\R) = I d'hd'hB^'\k,,k,) W{hR)WihR)W{\k, + k,\R). (45) 

For large scales, the expansion in Eq. (^) can be rewritten as = S^^ + S^^^a'^ + . . ., where: 

(R) = S^^^ (R) - 2 (R) (R) . (46) 
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The tree-level skewness has been thoroughly studied (Goroff et al. 1986, Juszkiewicz et al. 
1993, Bernardeau 1992, Bernardeau 1994, Lokas et al. 1995) including smoothing effects for 
both top-hat and Gaussian smoothing. The result for scale-free initial power spectra is: 

5f = y-(n + 3), (47) 
for top-hat smoothing (Bernardeau 1992b, Juszkiewicz et al. 1993), and: 



S. 



(0) 



3 



^ /n + S n + 3 3 1\ / 8\ ^ /n + 3 n + 3 5 1\ , 

''M— — ri)-{"+r)M— — -r4)' 



for Gaussian smoothing (Lokas et al. 1995, Matsubara 1994), where 2F1 denotes a hy- 
pergeometric function. More specifically, S^^ = 4.02, 3.71, 3.47, 3.28, 3.14 for n = 
— 2, — 1.5, — 1, — 0.5, for Gaussian smoothing. See Appendix P for the corresponding re- 
sults in the Zel'dovich approximation. 



3.2. Self-Similarity and Perturbation Theory 



Since there is no preferred scale in the dynamics of a self-gravitating pressureless perfect 
fluid in an Einstein-de Sitter universe, Eqs. (|l|)-(0) admit self-similar solutions (Peebles 
1980). This means that correlation functions of the cosmological fields should scale with a 
self-similarity variable, given appropriate initial conditions: knowing a statistical quantity at 
a given time completely specifies its evolution. For Gaussian initial conditions and scale-free 
power spectra, one can define a physical scale Rq, the correlation length (see Eq. (PP])), which 
obeys Rq linear theory (and in general if the non-linear power spectrum evolves 

self-similarly) . Statistical quantities in linear perturbation theory evolve self-similarly with 
i?o, e.g.. 



R-'P^'\k,T)^V'~'\kRo), (49) 

Vi?W(ki,k2,r) = B^^\k,Ro,k,Ro). (50) 

When loop corrections are taken into account, however, self-similarity may be broken by 
the appearance of new scales required by infrared and ultraviolet divergences in the loop 
integrations. In fact, one may consider a linear power spectrum Pii{k, r) given by a truncated 
power-law. 
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[ otherwise, 

where A is a normahzation constant; the infrared and ultraviolet cutoffs e and are imposed 
in order to regularize the loop integrations (SFl). In a cosmological N-body simulation, 
they would correspond roughly to the inverse comoving box size and lattice spacing (or 
interparticle separation) respectively. 

In the absence of the cutoffs, the spectrum (^) would be scale-free. The introduction 



of fixed (time- independent) cutoff scales e and kc in the linear power spectrum (^) breaks 
self- similarity, because they do not scale with the self-similarity variable kRo. The extent 
to which one can take the limits e ^ and kc ^ oo will determine whether one recovers 
self-similar scaling for the statistical properties of the density field. Infrared divergences, 
regulated by e, arise in individual diagrams when n < — 1, due to the divergence of the rms 
velocity field at large scales (Jain & Bertschinger 1996, SFl). These divergences are just a 
kinematical effect and cancel when the sum over diagrams is done, as a consequence of the 
Galilean invariance of the equations of motion (SFl). 

Ultraviolet divergences, on the other hand, arise due to small-scale power, and become 
more severe as n increases. In fact, for n > —1, one-loop corrections to the power spectrum 
break self-similarity (SF2) 



(fci^o)" 61 {kR 



\2n+3 



27rr(^) 3157r(n+l)r2 



ra+3 
2 




where 77 = — (n + 1) is an exponent which measures the deviation from self-similar scaling, 
and the self-similarity breaking factor becomes a logarithm when n = —1. 

We have done a similar calculation for the "tadpole" diagram (-B||i and B^u) contri- 
butions to the bispectrum and found that the same self-similarity breaking factors appear 
in this case. For the other contributions, the explicit calculation is not possible to do ana- 
lytically, but it can be checked by numerical integration that the full one-loop bispectrum 
breaks self-similarity for n > —1. For equilateral configurations, this calculation can be 
summarized by the one-loop result for n > —1: 

B^'\kRo, kRo) = 6„ (fci?o)'"+' (^^J , (53) 
with bn a finite n-dependent constant factor and logarithmic terms breaking self-similarity 
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as n — > —1. This breaking of self-similar evolution, not seen in numerical simulations, is a 
feature of the perturbative approach; it comes from large loop momenta due to increasing 
small scale power as n increases. We therefore do not expect our results to be physical 
as n — > — 1 from below. On the other hand, for —3 < n < —1, one-loop corrections to the 
bispectrum scale as {kRo)'^'^^^ and preserve self-similar evolution. In this case, as was already 
considered for the power spectrum in SF2, it is more convenient to regularize the individual 
one-loop bispectrum diagrams by using dimensional regularization (see Appendix 0), which 
effectively takes the limits e — and A;c — * oo. We now turn to the results of this calculation. 



4. One-Loop Results: Entering the Non-Linear Regime 

4.1. Bispectrum 

We now consider one-loop corrections to the bispectrum for initial power-law spectra 
Pii{k) oc A;" with spectral index —3 < n < —1. In this case, the resulting bispectrum obeys 
self-similarity and, based on previous results for the power spectrum (SF2), the perturbative 
approach is expected to give a good description of the transition to the nonlinear regime. 
Due to statistical homogeneity and isotropy, the bispectrum -B(ki, k2, r) in the scaling regime 
(e -C fcj -C kc) only depends on time, the quantities ki, k2, and the angle 9 { ki ■ k2 = cos 6*). 
In order to display the analytic results, however, it is more convenient to trade the variable 
9 for the third side of the triangle, fcg = |ki ks]. Let B^^^ki, kg) = A^aV b^^\ki, k2, k^), 
with ki + k2 + ks = 0. Then, using the results of Appendix 0, and summing over diagrams 



according to the results in Section |3.1| , the one-loop correction to the bispectrum for n = —2 
reads: 



30279 2635 fci^ 37313 fci 38431 

i i 1 

34496 fci^ 51744/^2^ 206976 /ta^ 68992 A;i ^2^ 

233 fci^ 16517 197 A;/ 78691 fci^ 

^8624^7^57 ~ 362208^717 ^ 7392^7^7 ~ 275968^^^7 
23 A;/ 9791 fci^ 703 A;i^ 19867 fci^ 

i I i I t I i 

103488 A;2^ A;3^ 206976 A;2^ A;3^ 68992 A;2^ A;3^ 206976 A;2 fca^ 

5311 A;!^ 42983 A;i 131 A;i 28393 

H \ 1 1 

34496^2^3^ 362208 A;2 A;3^ 3696 A;2H3^ 19712 A;i A;2 A;3 

53973 A;i^ 108685 A;iA;2 59599 A:i^ 

^193T776'fc7A>7 ^ TsilCMA^ ^ 362208A7a7 

-(-permutations. (54) 
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A well-known approximation scheme that provides insight into the physics of the non-linear 
regime is the Zel'dovich approximation (ZA) (Zel'dovich 1970). We therefore also consider 
the perturbative expansion for the Z A dynamics (see Appendix ^) , and calculate one-loop 
corrections to the bispectrum and skewness factor as well (loop corrections to the bispectrum 
in the ZA have been recently considered for spectra with small scale cutoffs and spectral 
indices n > —1 by Bharadwaj (1996)). For the one-loop bispectrum, again for n = —2, we 
find: 



.^7 



di)n , , ^ 27 9h 27 2187fc 

^ ' ' ^' 128 256 A;2^ 256 fci A;2 32768 fca^ ^3 



9ki^ 2853 /ti^ 9ki^ 1467 fci^ 



+ ^w^^^^;;^! ^^^^^ 

256fc2^fc3^ 16384 A;2^ fcs^ 256 A;2^ fcg^ 32768/^2^3^ 
2493A;iA;2 3 A;/ 3ki'^ 4863 fci^ 

8192 A;3^ 256A;2^fc3^ 256 A;2 A;3^ 16384 A;2^ /^s^ 
3ki^ 237 ki 11823 

+ - ^I^^W + 16384 A;, A:2 ^3 + P^^^^^^^^' 

(55) 

Using the one-loop power spectrum for n = —2 given in SF2 (see also Makino et al. (1992)) 

55 

= 53^. (56) 

where P^^\k) = A^a'^n^ p^^\k), we can obtain the one-loop hierarchical amplitude Q^^^ 
from Eq. (pl). Since Q^^^ depends on time, a convenient parametrization of the degree of 
nonlinear evolution which takes advantage of self-similarity is to write wave- vectors in terms 



of the correlation length Rq defined in Eq. (29), which for scale-free initial power spectra 
and Gaussian smoothing gives 

R-^'^2nAa^T{^^). (57) 

Figures |^ and |] show the resulting hierarchical amplitude Q (see Eq. (|37D ), for the exact 
dynamics (ED) and the Zel'dovich approximation (ZA) respectively. We see that for the 
ED one-loop corrections to Q are in general not negligible even for weakly nonlinear scales. 
When kiRo ^ 1, the contribution to the variance per logarithmic interval A(A;) = A7ck^P{k) 
becomes of order one, and we expect one-loop perturbation theory to break down, since the 
scales considered become comparable to the correlation length. It is interesting that at these 
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Fig. 5. — The hierarchical amphtude Q (see Eq. (0)) for triangle configurations {ki = 1, 
k2 = 0.5; ki-k2 = cos 9) as a function of the angle 9 to one-loop. The lowest full curve shows 
Q at tree-level, whereas the subsequent curves correspond to different stages of non-linear 
evolution parameterized by the first side of the triangle in terms of the correlation length, 
kiRo = 0.1, 0.2, 0.4, 0.8, 1 (from bottom to top). 

scales Eq. (|37D saturates, that is, the one-loop quantities B'^^^ and S*^^) dominate over the 
corresponding tree- level values and further time evolution does not change the amplitude Q, 
because 5(1) and have the same scale and, by self-similarity, time-dependence. Note that 
for this initial spectrum, the one-loop correction enhances the configuration dependence of 
the tree-level bispectrum as the system evolves to the non-linear regime. This enhancement 
is stronger for the ZA, which is understandable in view of the tendency of this dynamics 
to produce highly anisotropic two-dimensional structures (pancakes). Note that the ZA 
underestimates the one-loop correction, in correspondence with the unsmoothed skewness 
(SFl) and power spectrum results (SF2). 

Based on results from N-body simulations, it has been pointed out by Fry et al. (1993) 
(see also Fry et al. (1995)) that for n = —1 nonlinear evolution tends to "wash out" 
the configuration dependence of the bispectrum present at the largest scales (and given by 
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Fig. 6. — Same as Figure |^ for the Zel'dovich approximation. 

tree-level perturbation theory), giving rise to the so-called hierarchical form Q ^ const in 
the strongly non-linear regime. One-loop perturbation theory must predict this feature in 
order to be a good description of the transition to the nonlinear regime. To study this, we 
integrated numerically the one loop bispectrum for different spectral indices to understand 
the transition from the behavior at n = —2 to the n = — 1 spectrum (for n ^ —2 the one- loop 
bispectrum can be represented in terms of hypergeometric functions of two variables (see 
Appendix |^)). Figure ^ shows the result of such a calculation for the exact dynamics, in 
terms of the one- loop hierarchical amplitude Q*-^^ (see Eq. (^1])) for spectral indices running 
from —1.6 to —1.3. Clearly one-loop perturbation theory predicts a change in behavior of 
the nonlinear evolution: for n < —1.4 the one-loop corrections enhance the configuration 
dependence of the bispectrum, whereas for n > —1.4, they tend to cancel it, in qualitative 
agreement with numerical simulations. Note that this "critical index" n^ ^ —1.4 is the 
same spectral index at which one-loop corrections to the power spectrum vanish, marking 
the transition between faster and slower than linear growth of the variance of density field 
fluctuations (SF2) (see also Makino et al (1992), Lokas et al. (1995b), Bagla & Padmanabhan 
(1996)). Figure^shows that the same situation arises in the Zel'dovich approximation, which 
has ric ~ —1.75. Figure ^ displays the one-loop correction to the power spectrum in both 
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Fig. 7. — The one- loop hierarchical amplitude Q^^^ (see Eq. (|41|)) for triangle configurations 
(fci = 1, /c2 = 0.5; ki ■ k2 = cos9) as a function of the angle 9 for different spectral indices 
n. The spectral index runs from n = —1.6 (top full curve) to n = —1.3 (bottom full 
curve) in steps of An = 0.02. This shows the transition from positive to negative one- 
loop corrections as n is increased. The linear power spectrum in this figure is such that 
Vn{k) = {kRorexp[-{e/ky]exp[~{k/k,y]/[27rT[{n + 3)/2]]. 



the exact dynamics (ED) and the Zel'dovich approximation (ZA) in terms of the function 
a(n) defined by 



r{kRo) = 



{kRo 



27rr[(n + 3)/2] 



1 + a(n) {kRo) 



n+3 



(58) 



obtained by dimensional regular izat ion in SF2. Figures |, and ^ clearly illustrate the 
special character of the critical index in both dynamics. Note that in Figures ^ and the 
calculation is done by numerical integration and the linear power spectrum is not exactly 
scale-free, which can account for the very small shift in the critical index in these figures 
with respect to the exact scale-free case in Fig. p|. 



The change in behavior of the one-loop corrections at n ~ can be understood in phys- 
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Fig. 8. — Same as Figure |^ for the Zel'dovich approximation. 

ical terms as follows. As we increase n, the increase in small-scale power generates random 
motions which tend to prevent the collapse of high density regions (the "previrialization" 
effect (Davis & Peebles 1977, Evrard & Crone 1992, Lokas et al. 1995b, Peebles 1990). This 
is reflected in the sign change of one-loop corrections to the variance, which measures the 
growth of fluctuations. Another manifestation of this effect, which influences the shape of 
the bispectrum, is that random motions due to small-scale power cause structures to be 
less anisotropic and flows to have a smaller coherence length. This leads, upon ensemble 
averaging, to a cancellation of the configuration dependence in the hierarchical amplitude 
Q in Fourier space. A quantitative comparison of the predictions of one-loop perturbation 
theory with N-body simulations for the hierarchical amplitude Q is under way and is the 
subject of a forthcoming paper (Scoccimarro et al. 1996). 
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Fig. 9. — One- loop corrections to the power spectrum in the exact dynamics (ED) and 
Zel'dovich approximation (ZA) as a function of spectral index (see Eq. (^)). 
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4.2. Skewness: Comparison with Numerical Simulations 

We now consider one-loop corrections to the skewness factor. Since tree-level one-point 
cumulants such as the skewness and kurtosis are given by the spherical collapse dynamics 
(Bernardeau 1992, Bernardeau 1994b), one-loop contributions contain the first corrections 
to the spherical model coming from tidal motions. Given the analytic results for the n = —2 
bispectrum in the previous section, we can use Eq. (|45| ) to integrate numerically the one-loop 
skewness for different window functions. For the exact dynamics, we obtain (see Eq. (^)): 



Ri 3.86 + 3.18 (59b) 



for top-hat smoothing, and 



for Gaussian smoothing. For the Zel'dovich approximation we get: 

for top-hat smoothing, and 

S3-(ii) = ^H±^^^« 3.14+ 1.00 (62) 

^ ^ ^ 1 + 0.59 (T^(i?) ^' ^ ^ 

for Gaussian smoothing. We also computed one-loop corrections to the skewness in the Local 
Lagrangian approximation scheme discussed by Protogeros & Scherrer (1996). We obtain 
(see Appendix |Cp : 



for top-hat smoothing. Note that the results in Eqs. (p9D-( ^72D are all for n = —2. A visual 
summary of the exact perturbative results is given in Fig. 0, where we compare to the 
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numerical simulations by Colombi et al (1996). These N-body simulations used a tree code 
(Hernquist, Boucliet & Suto 1991), with 64^ particles in a cubic box with periodic boundary 
conditions. Symbols in Fig. |TU| correspond to different output times: diamonds (a = 2), 
triangles (a = 3.2), and stars (a = 5.2). The estimated systematic uncertainties in these 
measurements of skewness is ibO.l in logarithmic scale, which is due to the uncertainty in 
the finite volume correction applied to the simulation data (see Colombi et al. (1996) and 
also Colombi et al. (1994), Hivon et al. (1995) for details). This correction, which is rather 
large because of the large-scale power present in the n = —2 spectrum, is more important 
when the correlation length becomes a non- negligible fraction of the box size (i.e., for larger 
a). For a given time output, the finite volume correction is more important for large scales. 
Note that only the two latest outputs in this figure have been corrected for finite volume 
effects; for the first output (a = 2) this correction should be negligible. 

The two solid curves in Fig. ^ correspond to the predictions given in Eq. (|59a|) (bottom) 
and Eq. ( ^9b|) (top). The lower curve shows a saturation at values of cr^ ~ 1 such that one- 
loop corrections in Eq. (|59a|) dominate over the tree-level contributions, similar to what 
happens with the hierarchical amplitude Q. This saturation value, however, is not in good 
agreement with the numerical simulation data, which is not surprising given that those scales 
are well into the non-linear regime. We note that the N-body results are systematically lower 
(although within the error bars) than the one-loop perturbative calculation in the weakly 
non-linear regime. In fact, as cr^ they approach asymptotically the tree- level value given 
in the Zel'dovich approximation. This is most likely an artifact coming from the fact that 
the simulation uses ZA initial conditions, which have not been erased by the relatively early 
output time at a = 2 (Baugh et al. 1995). Note that the dashed curves given by: 



9.97 a^H{R) 



1 + 1.76 (T^niR)' 
3 + 4.69 4^(i?), 



(64a) 
(64b) 



which are the tree-level value given by the ZA plus the one-loop correction in the exact 
dynamics, fit the numerical results better, suggesting that indeed transients from the ZA 
initial conditions are still present in the first output. It is interesting to note that the 
expansion for large scales given in Eq. (|64b| ) seems to describe the transition to the non- 
linear regime better than Eq. ( |64a| ) , which soon becomes dominated by one-loop corrections 
and driven to the saturation value. Overall, however, we see that one-loop perturbation 
theory agrees with the simulation within the error bars even on scales where o"^ ~ 1, and 
therefore describes most of the transition from the tree-level result (valid in the limit cr^ 0) 
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Fig. 10. — The skewness factor as a function of the variance of density fluctuations cr^(i?) 
for spectral index n = —2. The symbols show the results from numerical simulations by 
Colombi et al. (1996) for top-hat smoothing. Different symbols correspond to different 
output times: diamonds (a = 2), triangles (a = 3.2), and stars (a = 5.2). Estimated error 
bars in these measurements are ±0.1 (systematic) in logarithmic scale (Colombi et al. 1996). 
The solid curves correspond to the prediction for top-hat smoothing (TH) of exact one-loop 
perturbation theory (ED), Eqs. ( |59a|) (bottom) and ( ^9b|) (top). Dashed lines denote the 
tree-level Zel'dovich approximation (ZA) plus one-loop ED, Eqs. ( |64a| ) (bottom) and (|64 
(top). Dot-dashed lines correspond to tree- level values in ED and ZA. 



to the nonlinear regime where 5*3 approximately approaches a constant value, in agreement 
with the corresponding results for the power spectrum (SF2). A more detailed comparison, 
with more accurate N-body measurements, will be presented elsewhere. 

The Zel'dovich approximation for 5*3 clearly underestimates the numerical simulation 
and the exact dynamics perturbative results, in agreement with previous results for un- 
smoothed flelds in SFl. Note that it also fails to describe properly the transition to the 
non-linear regime. The Local Lagrangian approximation (LLA) does reasonably well, al- 
though it overestimates the exact perturbative results; this is not surprising, given that it 
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has been "designed" to reproduce the tree-level one-point cumulants (Protogeros & Scher- 
rer 1996). On the other hand, a phenomenological model recently proposed by Munshi & 
Padmanabhan (1996), which assumes the hierarchical ansatz (equivalent to a constant 5*3), 
predicts S3 = 11.76 for —0.6 < log^o ^ 0.3, in disagreement with the numerical simulation 
results, which do not support the hierarchical assumption in the transition to the non-linear 
regime. 

It is interesting to note the importance of smoothing in the value of one-loop corrections 
by comparing the above results to the unsmoothed values found in SFl. For the exact 
dynamics, we obtained for unsmoothed fields 



^ ^ 4.86 + 10.03 a', (65) 

for the n = —2 spectrum, where the corresponding result for the Zel'dovich approximation 
reads: 



S. 



ZA 



4 + 4.69 



(66) 



In each case, smoothing reduces the relative importance of the loop corrections. For the 
exact dynamics, one-loop corrections to the smoothed skewness begin to dominate over the 
tree- level contribution for o"^ ~ 1, instead at cr^ ^ 1/2 for unsmoothed fields. 

Another interesting issue is the spectral dependence of the one- loop corrections to S3. 
One-loop corrections to the variance and average two-point correlation function show a linear 
dependence on spectral index for — 3 < n < — 1.5 (SF2), and we conjecture that a similar 
behavior extends to 5*3. For n = —3, smoothed and unsmoothed quantities coincide, because 
small-scale filtering does not affect statistical properties for a model with such extreme large- 
scale power (e.g., the variance is infrared-divergent). In this case, we have (SFl) 



Ss{R) 



34 

y 



10 (T^HiR) 



[n 



-3). 



(67) 



Taking into account the above results for n = 
with n, we expect that for top-hat smoothing. 



-2, Eq. 



9|) and assuming linear behavior 



■34 , 
- + 10 4^ 



a^HiR) 10 - 6.8 {n + 3) 



(R) 



— (n 



3) 1 + 6.8 (x^HiR) 



(68) 
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In fact, as n ^ —1.6, this ansatz leads to 5*3^^ — > 0, which is probably outside the region of 
validity of the linear extrapolation. For n — —1, for example, numerical simulations show 
that > (Colombi et al. 1996). 



5. Conclusions 

We have calculated one-loop corrections (non- linear corrections beyond leading order) to 
the bispectrum and skewness of the cosmological density field including smoothing effects, 
induced by gravitational evolution for Gaussian initial conditions and scale-free initial power 
spectra. These results extend previous calculations done at tree-level and allow us to probe 
the transition to the non-linear regime. 

We have shown that the one-loop bispectrum follows a similar behavior as a function of 
spectral index as the one-loop power spectrum. For n < — 1.4, one-loop corrections increase 
the configuration dependence of the bispectrum; for n > — 1.4, one-loop corrections tend to 
cancel the configuration dependence of the tree-level bispectrum, in agreement with n = —1 
numerical simulations. Therefore, there is a "critical index" ^ —1.4, similar to what 
happens in the case of the power spectrum, where one-loop corrections become negative 
at n > —1.4 (SF2), indicating a slowing down of the growth of fluctuations (Lokas et al. 
1995b). This increase in the configuration dependence of the bispectrum for n < is a 
prediction of one-loop perturbation theory that can be tested against numerical simulations. 
The configuration dependence of the bispectrum is due to the anisotropy of structures and 
flows in real space, and therefore has a direct physical meaning. In this respect, the one- 
loop bispectrum for the Zel'dovich approximation, which is well known to produce highly 
anisotropic structures (pancakes), shows a stronger configuration dependence, as expected. 
We interpret the change in behavior of the bispectrum as the spectral index increases as a 
result of the increased effect of small-scale power in the collapse of high density regions (the 
"previriahzation" effect (Davis & Peebles 1977, Evrard & Crone 1992, Lokas et al. 1995b, 
Peebles 1990). The random motions due to small-scale power tend to slow down the collapse 
and disrupt coherent structures and flows on small scales, which is reflected in the one-loop 
corrections to the power spectrum and bispectrum respectively. 

For spectral indices n < —1, self-similarity is retained at the one- loop level, and one 
can calculate loop corrections in the scaling regime by using the technique of dimensional 
regularization. We obtained explicit analytic results for the one-loop bispectrum for n — —2 
initial conditions; this allowed us to calculate the skewness factor for top-hat and Gaussian 
smoothing. We then extended the results in SFl to include local averaging of the fields; for 
n = —2, the smoothed one- loop correction is reduced by more than a factor of 2 from its 
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unsmoothed value, which shows the importance of smoothing in determining the value of one- 
loop corrections. The results for top-hat smoothing compare well with the corresponding 
measurements in numerical simulations, providing a description of most of the transition 
from the tree-level value at large scales {a 0) to the non-linear regime where S3 attains 
an approximate constant value in accord with expectations based on stable clustering. 

The results presented in this work suggest future directions in which one can improve 
the current understanding of the transition to the non-linear regime. One obvious extension 
would be to consider spectral indices n > —1, to generalize the present results for arbi- 
trary scale-free spectra. This would involve taking into account the effects of small-scale 
fluctuations on the evolution of large-scale modes in a way that absorbs the divergences 
that appear in the present formalism (renormalization) , and therefore recovers self-similar 
evolution for statistical quantities such as the power spectrum and bispectrum. This would 
lead to a better understanding of the role of previrialization in determining the structure of 
the correlation functions on intermediate scales. We hope to come back to this point in the 
near future. 

Since realistic power spectra are not scale-free, an important further step is to consider 
initial conditions such as those given by the cold dark matter (CDM) model and its variants. 
Since these models have effective spectral indices in the range rie// ~ — 2 to — 1 over the 
scales of interest, we expect that they will show similar features to the ones we presented 
here. Nevertheless explicit calculations are required in order to assess the effect of one-loop 
corrections in the determination of bias from the configuration dependence of the bispectrum 
(Fry 1994). Similarly, recent claims by Jing & Borner (1996) that there is a discrepancy in 
the three-point function between tree-level perturbation theory and numerical simulations 
for CDM models, may be properly addressed by taking into account one-loop corrections. 
As we showed in this work, these can be non-negligible even on weakly non-linear scales, 
depending on the initial spectrum. Work is in progress on these issues (Scoccimarro et al. 
1996). 

There is clearly much more work to do to understand non-linear clustering in an expand- 
ing Universe. However, the interplay between perturbation theory and N-body simulations 
suggests that there are three distinct regimes that describe its most important statistical 
features. At the largest scales, tree-level perturbation theory is well established as providing 
a good description of the correlation functions and the a — limit of the Sp parameters (Fry 
1984, Bernardeau 1994). In the strongly non-linear regime, numerical simulations (Hamil- 
ton et al 1991, Peacock & Dodds 1994, Suto 1993, Jain et al. 1995, Colombi et al. 1996, 
Jain 1996) have shown reasonable agreement with the stable clustering hypothesis, although 
there are still large uncertainties due to limitations in dynamic range. In this regime, there 
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is as yet no compelling analytic model which makes predictions in good agreement with the 
numerical results, except probably for the two-point function (Shcth & Jain 1996). In par- 
ticular, there is no understanding of the hierarchical structure of the Sp parameters, which 
seem to reach a plateau in the highly non-linear regime, Sp ~ constant. Finally, the results 
presented in this work suggest that the third regime, the transition to the non-linear regime, 
with cr ~ 1, can be understood by one-loop perturbation theory for models without excessive 
small-scale power. This is clearly promising and deserves further investigation. 
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S. Meyer, and M. Turner for comments and discTissions. Financial support by the DOE at 
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A. Dimensional Regularization 

To obtain the behavior of the one- loop N-point spectra forn < —1, one can use dimen- 
sional regularization (see e.g. Collins (1984)) to simplify considerably the calculations. Since 
we are interested in the limit kc ^ oo, all the integrals run from to oo, and divergences are 
regulated by changing the dimensionality d of space: we set d = 3 + e and expand in £ -C 1. 
For the bispectrum, we need the following one-loop three-point integral: 

^('"■''-''')"/(4=M(k,-q?M(k.-q)T3- 

When one of the indices vanishes, e.g. = 0, this reduces to the standard formula for 
dimensional- regularized two-point integrals (Smirnov 1991): 

r(i/i)r(i/2)r(d- t/i - 1/2) 

The integral J(i^i, 1^2, i^s) appears in triangle diagrams for massless particles in quantum field 
theory, and can be evaluated for arbitrary values of its parameters in terms of hypergeometric 
functions of two variables (Davydychev 1992). The result is: 
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r(z/3)r(z/i23 - d/2) 



r(z/i)r(z/2)r(z/3)r(d - uus) \ 

xF4(z/3, z/123 - d/2; 1 + U23 - d/2, 1 + z/ia 
xT{d/2-urs)T{d/2-U2s)+y'^'-''''T{u._^ _ 
xF4(z/2, d/2 - z/i; 1 + z/23 - c//2, 1 - U13 + d/2; x, y) 
xr(z/i3 - d/2)T{d/2 - z/23) + x'^/2-^^^r(z/i)r(rf/2 - z/2 
xF4(z/i, d/2 - z/2; 1 - «^23 + c?/2, 1 + z^i3 - d/2; x, y) 

xr(rf/2 - Z/i3)r(z/23 - rf/2) + X°'/2-<^23^«!/2-mr^^/2 

^i{d - z/123, rf/2 - z/3; 1 - z/23 + d/2, 1 - z/13 + d/2; x, y) 

v{d - z/i23)r(z/23 - f//2)r(z/i3 - rf/2)^ , 



(A3) 



where z/i23 = z/i + z/2 + z/3, Vij = z/j + z/^, a; = (k2 — ki)^/fc^, y = k^/kl, and F4 is Apell's 
hypergeometric function of two variables, with the series expansion: 



p . , , ^ ^^^'y' («).+i(&).+. 
F4(a,6; c,d;x,y) = 2^}^ 



i=0 j=0 



(c)i(rf)j 



(A4) 



where (a)j = T{a + i)/T{a) denotes the Pochhammer symboL When the spectral index is 
n = —2, the hypergeometric functions reduce to polynomials in their variables due to the 
following useful property for —a a positive integer: 



F4 a, b; c,d;x,y) = -— J -. 

1^0 J^o j! {c)i{d)j {-a-i-j)\ 



(A5) 



When using Eq. (|A3|), divergences appear as poles in the gamma functions; these can 
be handled by the following expansion (n = 0, 1, 2, . . . and e —* 0): 



T(-n + e) 



^ + 7/>(n + 1) + I + 7/.2(n + 1) - ij\n + 1) j + 0{e^) 



(A6) 



where ip{x) = d\iaT{x)/dx and 



1 1 

ij{n + l) = 1 + - + ... + - -7^ 

2 n 



(A7) 
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k=i 



with ipil) = -7e = -0.577216 ... and V^'(l) = 71^6. 



B. Zel'dovich Approximation 

In this approximation (Zel'dovich 1970, Shandarin & Zel'dovich 1989), the motion of 
each particle is given by its initial Lagrangian displacement. In Eulerian space, this is 
equivalent to replacing the Poisson equation by the ansatz (Munshi & Starobinski 1994): 

v(x,r) = -^^V$(x,r), (Bl) 

which is the relation between velocity and gravitational potential valid in linear theory. The 
important point about the ZA is that a small perturbation in Lagrangian fluid element paths 
carries a large amount of non-linear information about the corresponding Eulerian quantities, 
since the Lagrangian picture is intrinsically non-linear in the density field. This leads to non- 
zero Eulerian perturbation theory kernels at every order. The ZA works reasonably well as 
long as streamlines of fiows do not cross each other. However, multistreaming develops at 
the location of pancakes, leading to the breakdown of ZA (Shandarin & Zel'dovich 1989). 
The equations of motion in Fourier space are: 



d6{k, r) 
dr 



e{k, r) = -J (fki J d^k26D{k - ki - k2)a(k, ki)^~(ki, r)5(k2, r), (B2) 



d9(k,T) n{T) ~ 



dr 



e{k,r) = -Jd^k,Jd^k26D{k-ki-k2)P{kMM)0{ki,r)e{k2,r) . (B3) 



These equations, together with the perturbative expansion (|^), lead to the recursion relations 
in > 2): 



n-l 



-^jf (qi) • • • ) Qn) — X! • • • ' ^" 



m=l 



«(k,ki)^ . X /5(k,ki,k2) 



n 



n[n 



XG„_^(qm+l, • • • , ^.n] 



(B4) 
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Gni<il, • • • ,qn) = X! ^m(^l' • • • ' ^m) / ' _\\ Cn-mlcim+l) • • • y ^n) , (B5) 

m=l l'^ 

After symmetrization, the density field kernel takes the simple form (Grinstein & Wise 1987): 

n\ qf ql 

Using Eq. ( |B6| ) for n = 2 and (0), one can can calculate the skewness at tree- level for 
different window functions. For top-hat smoothing, the result is (Bernardeau 1994) 

5f=4-(n + 3), (B7) 

whereas for Gaussian smoothing, using the methods described by Lokas et al. (1995), we 
obtain: 



2 ' 2 ' 2' 4^ V 7 z 2 ' 2 ' 2' 4 
{n + Sy /n + 5 n + 5 7 1 



30 2 H 2 ' 2 '2'4 

where 2F1 denotes a hypergeometric function. More explicitly, for Gaussian smoothing, 
= 3.14, 2.80, 2.51, 2.26, 2.04 for n = -2, -1.5, -1, -0.5, respectively. 



Fi— ^,-,T, B8 



C. Local Lagrangian Approximation 

In this approximation (Protogeros & Scherrer 1996), the final density at a Lagrangian 
point q at time r is assumed to be a function only of the initial density at the same Lagrangian 
point and time r: 

/ N_ ^o(q) 

"''>'^)=ll-oM«„(q)/„l- 

where 1] = 1+5, with So{q) = 5(q, tq) and tq is the initial time. The constant a takes the value 
a = 1 for the planar approximation (which becomes exact for one-dimensional collapse), 
whereas a = 3 corresponds to spherical collapse in the Zel'dovich approximation. Local 
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Lagrangian approximations at the level of the equations of motion for fluid elements have 
been considered recently by Hui & Bertschinger (1996). In this work we restrict ourselves 
to the case a = 3/2, which, although it has no particular physical meaning, can be shown 
to closely approximate the hierarchical amplitudes Sp of tree-level perturbation theory for 
the exact dynamics (Bernardeau 1992). Upon normalization of the probability distribution 
function for 5, the local Lagrangian approximation in this case reads: 

([l-2 5i(q,r)/3p/2)^ 

[l-2^,(q,r)/3]3/^ ' (^'^ 

where <>iag denotes ensemble averaging in Lagrangian space, and 5i(q, r) corresponds to 
the evolution of the density contrast in linear perturbation theory. The 5*^ parameters are 
defined as > 2) 

= < §2 >p-i ' (C3) 

where the angular brackets correspond to Eulerian ensemble averages. For powers of rj we 
have (Protogeros & Scherrer 1996) 



< ry- >=< r^rn-i _ (C4) 

Given that, for Gaussian initial conditions, < 5™ >= ("^ ^ 1)- '^^ ^ even and zero 
otherwise, one has everything needed to compute Eq. (|G3[ ) for unsmoothed density fields. 
The effects of top-hat smoothing for power-law power spectra can be included by the implicit 
mapping (Bernardeau 1994b, Protogeros & Scherrer 1996) 

where /(x) = (1 — 2x/3)^^^'^ denotes the unsmoothed mapping. We are particularly inter- 
ested in n = —2, for which Eq. ( |C5| ) yields: 

r]s = z-\6i) < z%6i) >iag , (C6) 

where ^(^i) denotes the appropriate solution to the quartic equation = 1 — 2Siz/3. Using 
Eqs. ( |C3D and (|U4D , we obtain for the skewness and kurtosis factors 
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4 + 509/54 4^^4 ,209 2 -4 + 387 n-^ (R) (Cr) 
^'^^^ = 1 + 25/18 aUR) ^ ^M- ^^^^^^ ^ ^ + ^-^^ ^^^^ 



269/9 + 53661191/373248 269 30419591 ^ 

~ 1 + 25/12 4^ (i?) ^"9~+ 373248 

^ 29.88 + 81.49 4h(^)- (C^ 
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